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Abstract  The  stress-strain  behavior  of  cross-linked  polymeric  networks  was  investigated 
using  molecular  dynamics  simulations  with  a  coarse-grained  representation  of  the  repeating 
units.  The  network  structure  was  formed  by  dynamically  cross-linking  the  reactants  placed 
between  two  rigid  layers  comprised  of  particles  of  the  same  type.  We  studied  two  types 
of  networks  which  differ  only  by  one  containing  ionic  pairs  that  amount  to  7%  of  the  total 
number  of  bonds  present.  The  stress-strain  curves  were  obtained  after  imposing  deformation 
in  tensile  and  shear  modes  to  the  networks  and  measuring  their  stress  response.  Under  both 
forms  of  deformations  there  was  improvement  in  the  level  of  stress  that  the  material  could 
bear.  Moreover,  the  time  dependent  behavior  of  the  improvement  in  mechanical  properties 
signified  a  self-healing  mechanism. 

Keywords  Polymer  networks  •  Self-healing  •  Ionic  bonds  •  Toughness  •  Molecular 
dynamics  simulations  •  Coarse-grained  modeling 


1  Introduction 

One  of  the  challenges  that  material  science  is  facing  is  to  design  smart  materials  that  can 
sense  the  presence  of  a  defect  and  actively  repair  the  defected  area.  Such  self-healing  materi- 
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als  would  significantly  extend  the  lifetime  and  the  usefulness  of  manufactured  items  (Balazs 
2007)  such  as  polymer  matrix  composites  as  replacements  of  metals,  and  cross-linked  sys¬ 
tems  as  adhesives  or  coatings. 

The  approaches  addressing  the  self-healing  phenomenon  in  materials  science  can  be  cat¬ 
egorized  in  three  classes:  the  healing  with  an  embedded  liquid,  the  thermally  activated  solid 
phase  healing  and  the  healing  of  projectile  puncture  (Kessler  2007).  The  first  approach  in¬ 
volves  embedding  capsules  containing  a  reactive  healing  liquid  (e.g.,  resin)  that  rupture  and 
release  the  liquid  upon  deformation.  While  it  has  proven  to  be  successful  in  repairing  the 
damaged  area,  it  has  stability  issues  that  occur  during  processing  and  the  available  liquid 
payload  capacity  is  often  limited  for  various  reasons.  The  second  approach  provides  an  im¬ 
proved  performance  after  every  cycle  of  heating,  but  there  is  a  requirement  of  external  heat 
supply  that  makes  the  process  non-autonomic.  The  third  approach  has  gained  recent  inter¬ 
est  perhaps  because  it  neither  requires  an  external  intervention  nor  has  stability  issues.  It 
has  been  shown  that  the  target  film  (based  on  ionomers  such  as  copolymers  of  ethylene 
and  methacrylic  acid)  that  was  penetrated  by  a  projectile  with  a  high  speed  is  healed  al¬ 
most  instantly  to  an  airtight  condition.  Essentially,  the  penetration  of  the  projectile  causes 
localized  heating  near  the  puncture  which  then  closes  upon  itself  to  heal.  The  mechanism 
behind  such  a  self-healing  process  after  the  ballistic  penetration  involves  two  transforma¬ 
tions.  First,  the  film  which  has  been  locally  melted  due  to  the  heat  dissipation  re-cry stallizes 
and  solidifies.  Later,  the  reordering  of  the  physical  cross-links  takes  place.  These  physical 
cross-links,  according  to  Eisenberg-Hird-Moore  model  (Eisenberg  et  al.  1990),  are  primary 
aggregates  that  consist  of  several  ion  pairs  known  as  multiplets.  These  multiplets  bring  sev¬ 
eral  chains  together  bound  by  ionic  interaction  creating  physical  cross-links.  Calorimetric 
studies  (Tadano  et  al.  1989)  have  shown  that  the  re-crystallization  occurs  rapidly,  while  the 
reformation  of  the  physical  cross-links  is  a  rather  slow  process.  As  a  result,  the  material  fully 
recovers  after  the  latter  process  is  complete. 

Highly  cross-linked  epoxies  or  polyurethanes  are  a  class  of  materials  that  can  particularly 
benefit  from  a  self-healing  mechanism.  While  they  possess  a  variety  of  appealing  properties, 
two  major  shortcomings  of  the  highly  cross-linked  materials  are  their  brittleness  and  low 
elongation,  which  are  due  to  the  limited  degree  of  freedom  available  to  the  strands  between 
the  cross-linking  points.  The  strands  extend  to  their  limits  even  with  a  small  deformation 
and  cause  breaking  of  covalent  bonds.  Since  the  covalent  bonds  are  non-reformable,  any 
damage  that  is  caused  by  impact  or  elongation  will  be  permanent  and  cause  formation  of 
micro-cracks.  For  instance,  cyclic  mechanical  loads  result  in  accumulation  of  the  micro¬ 
cracks  which  then  leads  to  catastrophic  failure. 

The  important  question  is:  would  a  similar  idea  to  that  behind  the  healing  of  projectile 
puncture  in  the  thermoplastic  ionomers  apply  to  the  thermosets?  It  is  known  that  the  break¬ 
ing  of  covalent  bonds  is  responsible  for  reducing  the  impact  resistance  properties.  The  in¬ 
herent  toughness  of  the  network  is  accommodated  mainly  by  non-covalent  interactions  such 
as  the  van  der  Waals  and  the  hydrogen  bond  interactions.  Their  effectiveness  is  due  to  their 
capability  of  breaking  and  displacing  under  stress  to  absorb  impact  energy.  However,  under 
high  energy  impact  conditions  such  weak  interactions  may  not  provide  a  sufficient  level  of 
strength  in  the  material.  Therefore,  in  theory,  incorporation  of  stronger,  yet  reversible  ionic 
bonds  (i.e.  two  neighboring  oppositely  charged  atoms)  into  the  polymeric  matrix  could  be 
an  alternative  and  effective  approach  to  improve  the  impact  properties.  However,  one  needs 
to  realize  that  a  similar  mechanism  explained  above  for  the  thermoplastic  ionomers  is  an  un¬ 
likely  means  for  the  highly  cross-linked  materials.  First,  there  will  not  be  re-crystallization 
since  the  low  degree  of  freedom  in  the  material  prevents  crystallization.  Secondly,  the  for¬ 
mation  of  multiplets  is  hindered  also  due  to  the  cross-linked  nature  of  the  polymer. 
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For  the  cross-linked  systems  with  embedded  ionomers  one  could  envision  an  alternative 
mechanism  in  which  the  multiplets  are  replaced  by  ionic  pairs  and  the  extent  of  repairing 
action  is  reduced  primarily  to  the  length  scale  of  strands  between  cross-links.  A  temporary 
loss  of  the  ionic  interaction  caused  by  the  deformation  can  be  restored  by  the  same  or  another 
pair  of  ionic  groups.  Moreover,  the  addition  of  ionic  pairs  which  are  stronger  than  the  van 
der  Waals  and  the  hydrogen  bond  interactions,  yet  still  reversible,  could  serve  to  increase  the 
bulk  material  strength  and  the  toughness.  The  soundness  of  this  probable  mechanism  can 
suitably  be  studied  using  molecular  modeling.  Past  molecular  simulation  studies  of  ionic 
polymers  focused  mainly  on  their  charge  transport  properties  associated  with  applications 
such  as  poly  electrolyte  membranes  in  fuel  cells  (Din  and  Michaelides  1998;  Elliott  et  al. 
1999;  Vishnyakov  and  Neimark  2000).  Accordingly,  these  modeling  studies  have  primarily 
concerned  with  the  solubility,  the  conductivity,  and  the  diffusivity  rather  than  mechanical 
properties. 

Computational  modeling  of  networked  polymers  has  gained  limited  attention  in  contrast 
to  their  widespread  use  in  many  applications.  The  atomistic  models  of  such  systems  include 
epoxy  resins  that  were  modeled  using  a  single  step  cross-linking  procedure  (Yarovsky  and 
Evans  2002)  and  poly(dimethylsiloxane)  (PDMS)  network  where  a  dynamical  cross-linking 
procedure  was  used  (Heine  et  al.  2004).  Besides  these  atomistic  models  Stevens  and  cowork¬ 
ers  modeled  polymeric  networks  using  a  coarse-grained  representation  with  the  dynamical 
cross-linking  procedure  (Stevens  2001;  Tsige  et  al.  2004;  Tsige  and  Stevens  2004).  They 
investigated  the  effects  of  the  cross-linker  functionality,  the  network  connectivity  on  the 
mechanical  properties  and  the  influence  of  bond  density  on  interfacial  fracture.  Their  work 
has  shown  that  a  coarse-grained  representation  enables  simulation  of  the  networked  sys¬ 
tems  at  longer  time  and  length  scales  while  capturing  the  experimentally  observed  essential 
mechanical  features. 

Inspired  by  the  success  of  using  ionic  interactions  in  the  healing  of  projectile  puncture, 
here  we  seek  a  parallel  means  of  adding  self-healing  character  to  cross-linked  polymers. 
Specifically,  we  investigate  the  mechanical  properties  in  a  polymer  network  with  ion  pairs 
embedded  throughout  using  a  generic  coarse-grained  representation  of  the  system.  The  ionic 
pairs  were  incorporated  in  an  ideal  manner  as  explained  later  in  detail.  We  explore  the  pos¬ 
sibility  of  enhancing  the  toughness  by  ionic  links  and  the  potential  underlying  mechanisms 
using  molecular  dynamics  simulations.  We  have  paid  attention  also  to  the  effect  of  the  rate 
of  deformation  applied  to  the  network. 

This  article  is  organized  as  follows.  Section  2  describes  the  details  of  the  model  and 
the  simulation  protocol.  The  stress- strain  behaviors  of  the  standard  and  the  ionic  networks 
upon  deformations  in  the  tensile  and  shear  modes  were  presented  in  Sect.  3.  The  possible 
mechanisms  behind  the  results  were  discussed  in  Sect.  4.  The  conclusions  of  this  work  are 
given  in  Sect.  5  followed  by  the  appropriate  acknowledgments. 


2  Model  and  simulation  protocol 

The  cross-linked  system  was  represented  using  the  bead-spring  model  based  on  the  works  of 
Stevens  and  coworkers  (Stevens  2001;  Tsige  et  al.  2004;  Tsige  and  Stevens  2004).  The  inter¬ 
actions  between  the  beads  were  described  using  6-12  Lennard- Jones  (LJ)  potential  (see  (1)) 
with  a  2.5J  cutoff 


(1) 
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where  mq  is  the  LJ  energy,  d  is  the  LJ  diameter  and  r  is  the  distance  between  two  interacting 
beads.  Both  uq  and  d  were  taken  as  unity. 

The  beads  are  bonded  via  a  quartic  bond  potential  (see  (2)) 


k(y  -  bi)(y  -  b2)y^  +  Uq,  r  <rc 
Uq,  r  >rc 


(2) 


where  y  =  r  —  Ar  that  shifts  the  quartic  center  from  the  origin  and  the  parameters  as  adopted 
directly  from  the  study  by  Stevens  et  al.  (Stevens  2001)  are  k  =  \434.3uo/d^,  b\  =  — 0.759J, 
b2  =  0.0,  Ar  =  rc  =  l.5d  and  Uq  =  61 .223uq.  The  bond  potential  is  smoothly  cutoff  at  r^. 
At  distances  larger  than  the  bond  potential  is  turned  off  that  leaves  only  LJ  term  remain¬ 
ing.  The  bond  breakage  is  irreversible.  The  bond  potential  parameters  result  in  maximum 
bonding  force  of  156. compared  to  maximum  LJ  force  of  2AuQld. 

Coulombic  interactions  for  pairs  of  ionic  beads  were  computed  with  a  cutoff  of  also 
2. 5 J (for  the  direct  sum) 

jBqiqj 


Ucoa\i^)  =  ksT  - 


(3) 


where  qt  is  the  charge  on  bead  /,  kg  is  the  Boltzmann  constant,  T  is  temperature,  s  is 
dielectric  constant  and  Ib  is  the  Bjerrum  length.  Ib  is  defined  as  the  separation  distance 
at  which  the  electrostatic  interaction  between  two  interacting  elementary  charges,  in  a 
dielectric  medium  with  dielectric  constant  s  is  equivalent  to  the  thermal  energy 


h 


sUbT 


(4) 


Ib  was  set  to  50 J.  With  this  value  of  the  ratio  of  the  maximum  quartic  bonding  force 
to  the  ionic  interaction  force  becomes  about  4.  Specifically,  at  r  =  \.\2d  the  ionic  force  is 
39.1wo/^  versus  156. 7mo/^  (i-e.  ratio  of  4.0).  This  ratio  is  simply  an  estimate  value  that 
relates  the  strength  of  a  typical  covalent  bond  to  an  ionic  bond. 

The  MD  simulations  were  performed  using  LAMMPS  (Large-scale  Atomic/Molecular 
Massively  Parallel  Simulator)  (Plimpton  1995).  The  equations  of  motion  were  integrated 
using  the  Verlet  algorithm  (Verlet  1967)  with  a  time  step  of  0.05  r,  where  r  is  the  LJ  time 
unit.  The  long-range  Coulombic  interactions  beyond  the  cutoff  (reciprocal  sum)  were  calcu¬ 
lated  using  the  particle-particle  particle-mesh  Ewald  (PPPM)  solver  with  a  precision  value 
of  1.0  X  10  Both  the  temperature  and  pressure  were  controlled  using  Nose-Hoover  for¬ 
malism  with  a  damping  constant  of  O.lr. 

The  initial  simulation  box  consisted  of  small  chains  with  two  and  three  beads  that  were 
sandwiched  between  two  rigid  layers  of  beads.  The  chains  with  three  beads  contained  one 
bead  designated  as  the  cross-linker  bead  for  the  subsequent  cross-linking  process.  The  layers 
were  made  of  (111)  plane  of  the  face-centered  cubic  (FCC)  lattice  structure  with  lattice 
spacing  of  1.2.  They  interact  with  all  beads  except  themselves.  The  system  consisted  of  a 
total  number  of  103,344  beads  of  chain,  cross-linker  and  layer  type,  where  the  ratio  of  the 
number  of  two-  and  three-bead  chains  was  determined  by  stoichiometry.  The  liquid  mixture 
was  first  equilibrated  in  canonical  ensemble  (NVT)  at  1.0  reduced  temperature  unit  (T*), 
followed  by  the  density  equilibration  in  NPT  (isobaric-isothermal)  ensemble.  The  pressure 
was  imposed  anisotropically  in  z-axis  while  no  pressure  scaling  was  employed  in  x  and  y 
directions. 

The  liquid  mixture  was  dynamically  cross-linked  after  the  equilibration.  Briefly,  a  vary¬ 
ing  distance  criterion  was  used  first  to  determine  the  candidate  pairs  (i.e.  a  bead  and  a  cross - 
linker  bead  with  a  functionality  of  four)  which  then  were  sequentially  cross-linked.  The 
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selected  separation  distance  ranged  from  \.\d  at  the  first  step  until  l.35d  at  the  final  cross- 
linking  step.  Each  step  consisted  of:  1.  the  determination  of  pairs,  2.  the  creation  of  bonds, 
and  3.  the  equilibration  of  the  system.  The  cross-linking  (first  to  the  beads  in  the  layers  and 
then  in  the  bulk)  was  carried  on  for  a  total  of  40  and  60  steps,  respectively.  The  number  of 
steps  for  the  first  cross-linking  stage  was  adjusted  to  create  adequate  cross-linking  density 
on  the  layers  so  that  cohesive  failure  was  promoted.  The  second  cross-linking  stage  was 
continued  until  a  sufficient  degree  of  cross-linking  was  achieved  in  the  bulk.  Near  the  end 
of  each  cross-linking  stage,  only  a  few  or  no  additional  cross-linking  occurred.  As  a  result, 
about  60,000  cross-linking  bonds  were  created  with  nearly  5,000  of  them  bonded  to  the 
beads  in  the  rigid  layers,  leading  to  98%  degree  of  cross-linking.  After  equilibration  of  the 
network  at  LOT*,  the  harmonic  bond  potential  used  hitherto  was  replaced  with  the  quartic 
potential.  Subsequently,  the  system  was  slowly  cooled  down  (at  a  rate  of  1.4  x  10“^  T*  per 
MD  step)  to  0.3 r*,  which  is  well  below  the  estimated  glass  transition  temperature. 

Ionic  bonds  were  incorporated  into  the  fully  cross-linked  network  subsequent  to  its  equi¬ 
libration.  The  bond  between  a  pair  of  randomly  selected  consecutive  beads  was  removed 
to  create  ionic  bonds.  In  the  case  of  an  ionic  network,  these  beads  were  assigned  opposite 
charges  to  create  ionic  bonds  while  for  a  standard  network  they  were  left  unchanged.  Using 
this  approach,  we  created  a  new  network  with  the  same  network  density  (82%)  as  non-ionic 
one  but  with  ionic  bonds  incorporated  within.  The  fact  that  we  used  pairs  of  beads  bonded  to 
each  other  to  create  ionic  pairs  is  actually  expected  to  reflect  a  molecular  arrangement  close 
to  an  actual  system.  Namely,  one  would  anticipate  in  the  actual  synthesis  of  a  polymeric 
network  that  the  monomers  with  opposite  charges  would  be  at  close  proximity  as  the  cross- 
linking  reactions  take  place.  However,  there  could  also  be  multiplets  where  more  than  a  pair 
of  charged  particles  assembles.  Therefore,  our  approach  in  creating  the  ionic  network  from 
a  network  that  had  no  ionic  interactions  during  cross-linking  would  yield  an  ideal  network 
where  there  are  no  ionic  aggregates,  only  pairs. 

Two  modes  of  deformation  (tensile  and  shear)  were  applied  in  a  quasi-equilibrium  fash¬ 
ion.  In  the  tensile  mode,  the  upper  layer  was  displaced  along  the  z -direction  with  a  selected 
speed  over  an  increment  of  0.5  J.  Subsequently,  the  simulation  was  run  at  the  new  position 
for  the  same  number  of  simulation  steps  needed  for  the  displacement.  The  data  for  the  calcu¬ 
lation  of  the  stress-strain  curves  were  collected  at  this  stage.  The  effective  deformation  rate, 
therefore,  is  half  of  the  selected  displacement  speed.  The  two  selected  speeds  0.5  x  10“^ 
and  0.5  x  \0~^d/T  resulted  in  strain  rates  of  0.6  x  10“^  and  0.6  x  10“"^  1/r,  respectively. 
In  the  case  of  the  shear  deformation  the  upper  layer  was  displaced  in  the  v -direction. 


3  Results 

We  begin  our  analyses  with  an  inspection  of  the  microstructure  of  the  ionic  network  using 
a  set  of  radial  distribution  functions  (g(r))  shown  in  Fig.  1.  g(r) chain-chain  exhibits  an  initial 
sharp  peak  located  at  0.95  and  a  broader  peak  at  l.Old  corresponding  to  the  bonded  and  the 
non-bonded  interactions  among  the  beads  of  the  chain  type,  respectively.  In  other  words, 
the  quartic  bond  potential  with  the  parameters  used  in  this  study  shows  a  minimum  in  the 
potential  energy  at  0.93d  that  agrees  with  the  position  of  the  initial  peak,  while  the  energy 
minimum  in  the  Lennard- Jones  potential  for  the  non-bonded  interactions  is  at  1.12 J  in  line 
with  the  second  peak.  The  remainder  of  the  g{r) chain- chain  has  shallower  and  dying  out  peaks 
that  is  reminiscent  of  liquids.  g(r)^iinker-xiinker  is  very  similar  to  g{r) chain-chain  suggesting 
that  the  cross-linker  beads  are  uniformly  distributed  in  the  network.  However,  the  absence 
of  a  sharp  peak  near  0.93d  signifies  that  cross-linkers  are  not  bonded  to  one  another.  The 
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Fig.  1  Radial  distribution 
functions  of  pairs  of 
(a)  chain-chain  and 
xlinker-xlinker  and  (b)  ion-ion 
and  layer-xlinker  type  beads 


sharp  spike-like  peaks  in  g{r)xiinker -layer  indicate  a  layer  by  layer  positioning  of  the  cross¬ 
linker  beads  near  the  layers.  The  first  peak  is  a  result  of  cross-linkers  directly  bonded  to 
the  layers,  while  the  second  (~1.6J)  and  third  (~2.0(i)  peaks  are  due  to  the  cross-linkers 
that  are  separated  from  layers  by  one  and  two  layers  of  beads,  respectively.  The  intensity  of 
these  peaks  diminishes  rapidly  with  r.  This  type  of  arrangement  shows  that  the  network  has 
a  fairly  ordered  structure  near  the  layers  (r  <  ~5  J)  that  becomes  amorphous -like  away  from 
the  layers.  The  relative  positioning  of  the  ionic  pairs  are  represented  by  g(r) ton-ion  -  There  is 
a  single  and  dominant  peak  at  0.97 J  pointing  to  the  fact  that  they  were  created  from  pairs 
of  beads  which  were  initially  bonded  (i.e.  pairs  of  a  cross-linker  and  a  chain  type  beads). 

After  this  brief  examination  of  the  micro  structure  we  will  now  investigate  the  mater¬ 
ial  responses  to  the  tensile  deformation  at  low  strains  in  order  assess  whether  the  current 
model  follows  the  universal  relationships  between  modulus,  yield  stress  and  temperature  for 
cross-linked  systems.  The  stress-strain  curves  for  a  fully  cross-linked  network  (i.e.  98%  de¬ 
gree  of  cross-linking)  at  temperatures  ranging  from  0.3  to  0.557*  were  plotted  in  Fig.  2(a). 
The  Young’s  modulus  (7),  the  yield  stress  {Oy)  and  the  strain  at  yield  (^y)  were  extracted 
from  these  curves  and  are  presented  in  Table  1.  We  first  tested  the  Kambour’s  relationship 
(Kambour  1983), 

a,  =  C{T,  -  T)  (5) 

which  has  been  shown  to  be  generally  well  fitted  for  thermosets  (Cook  et  al.  1998; 
Tcharkhtchi  et  al.  1996).  The  linear  relationship  between  Oy  and  Tg  —  T  holds  for  our  model 
as  demonstrated  in  Fig.  2(b)  with  a  fitted  line  to  the  data  points. 

Next,  we  investigate  how  the  yield  stress  relates  to  the  Young’s  modulus.  It  is  well  known 
that  the  ratio  of  cfy/E  is  not  constant,  but  a  power  law  relationship  proposed  by  Kitagawa 
(1977)  was  found  to  be  valid  for  various  thermosets  (Cook  et  al.  1998;  Yamini  and  Young 
1980): 

2“^  =  (MV  (6) 

Tay,o  \TeJ 

where  Eq  and  Oy^o  are  the  modulus  and  the  yield  stress  at  a  reference  temperature  of  Tq. 
The  current  model  exhibits  the  same  power  law  relationship  between  Oy  and  E  as  presented 
in  Fig.  2(c).  However,  the  exponent  n  is  calculated  as  1.13  compared  to  the  experimental 
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Fig.  2  (a)  Stress-strain  curves  for  a  fully  cross-linked  network  (98%  degree  of  cross-linking)  at  varying 
temperatures  (T  =  0.3  (circle),  0.35  (square),  0.4  (diamond),  0.45  (triangle  up),  0.5  (triangle  left),  0.55 
(triangle  down)  T*).  (b)  Yield  stress  versus  Tg  —  T  with  a  linear  fit.  Tg  was  taken  as  0.5  T*  (Stevens  2001). 
(c)  Yield  stress  and  modulus  data  according  to  the  equation  proposed  by  (Kitagawa  1977).  Tq  =  0.47* 


Table  1  Young  s  Modulus,  yield 
stress  and  yield  strain  for  a  fully 

T(T*) 

E(uo/d^) 

cry(uo/d^) 

£y 

cross-linked  network  (98% 
degree  of  cross-linking)  at 

0.3 

50.0 

3.8 

0.090 

various  temperatures 

0.35 

43.8 

3.4 

0.095 

0.4 

40.6 

3.1 

0.096 

0.45 

34.9 

2.9 

0.098 

0.5 

30.8 

2.6 

0.100 

0.55 

27.0 

2.4 

0.110 

result  of  1.9.  The  difference  might  be  resulting  from  the  fact  that  the  deformation  rates  in 
the  simulations  are  several  orders  of  magnitude  higher  than  the  experimental  range.  Since 
the  strain  rate  appears  to  have  a  stronger  effect  on  the  yield  stress  than  the  modulus,  the  high 
deformation  rates  occurring  in  the  simulations  might  be  reducing  the  slope.  Another  reason 
could  be  the  low  chain  stiffness  in  our  model  (i.e.  the  absence  of  bending  and  torsional 
potentials). 

After  comparing  some  mechanical  properties  of  the  simulations  to  empirical  models  at 
low  strains  for  the  fully  cross-linked  network,  we  turn  our  attention  to  mechanical  responses 
in  the  tensile  and  the  shear  deformation  for  both  the  control  and  the  ionic  networks  over 
the  entire  strain  range.  The  stress-strain  curves  obtained  after  the  tensile  deformation  are 
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Fig.  3  Number  of  bonds  (iV^) 
(upper  panel)  and  stress  (cr) 
versus  strain  (s)  upon  the  tensile 
deformation  at  strain  rate  of 
0.6  X  10“^  1/r  for  the  standard 
network  (line)  and  the  ionic 
network  (circles).  Dashed  lines 
represent  data  for  a  network  with 
higher  degree  of  ionic  content 


shown  in  Fig.  3  (lower  panel).  There  is  a  clear  peak  for  both  networks  corresponding  to  the 
yield  behavior  at  s  0.1.  The  yield  stress  originates  primarily  from  the  LJ  interactions. 
For  this  reason,  the  yield  peak  appears  to  be  very  distinct  due  to  the  fact  that  the  networks 
have  a  lower  degree  of  cross-linking  than  a  fully  cross-linked  system  such  as  in  Fig.  2.  The 
degree  of  cross-linking  does  not  have  direct  influence  on  the  yield  stress;  rather  Oy  is  affected 
indirectly  due  to  the  change  in  Tg  (Pascault  et  al.  2002).  However,  the  plateau  is  reduced 
significantly  with  lowering  the  degree  of  cross-linking.  Therefore,  Oy  becomes  larger  than 
the  rest  of  the  stress-strain  curve  since  the  degree  of  cross-linking  was  lowered  to  82%. 
Note  also  that  the  equilibrium  distance  of  the  LJ  potential  is  at  1.12J,  while  the  maximum 
force  occurs  at  \.25d.  Therefore,  it  takes  a  strain  of  about  0.1  to  bring  the  network  from  the 
equilibrium  LJ  separation  to  a  separation  of  the  maximum  resistance. 

In  the  standard  network,  subsequent  to  the  yield  stress  there  is  plateau  region  that  lasts 
until  about  s  =  1.0.  During  this  period  the  strands  connecting  two  cross-linkers  are  being 
extended  while  bonds  remain  mostly  intact  as  described  by  Stevens  (2001).  The  bonds  start 
breaking  appreciably  only  after  about  s  =  1.2  (see  upper  panel  in  Fig.  3)  and  the  network 
fails  ultimately  at  about  s  =  2.0.  The  ionic  network  exhibits  a  similar  behavior  to  the  stan¬ 
dard  network,  but  with  some  imperative  differences.  First,  after  the  yield  behavior  the  values 
of  stress  do  not  drop  as  extensively  as  the  standard  network.  The  stress-strain  curve  remains 
above  that  of  the  standard  network  indicating  that  the  addition  of  ionic  bonds  improves  the 
toughness  (i.e.  the  area  under  this  curve).  Another  key  result  is  that  the  ionic  network  re¬ 
quires  nearly  twice  as  many  bonds  to  be  broken  to  fail  compared  to  the  standard  network. 
We  also  observe  that  the  ultimate  failure  does  not  occur  abruptly  due  to  the  lower  degree  of 
cross-linking. 

We  ran  a  set  of  simulations  for  a  network  with  14%  ions  to  observe  how  a  larger  ionic 
content  and  a  lower  degree  of  cross-linking  affect  stress-strain  behavior.  Following  the  same 
approach  as  in  the  Sect.  2,  the  resulting  network  had  66%  degree  of  cross-linking.  The  stan¬ 
dard  network  with  14%  ionic  content  loses  its  strength  (as  compared  to  7%  ionic  content) 
due  to  the  fact  that  the  network  loses  its  structural  integrity  at  the  low  level  of  cross-linking. 
In  the  case  of  the  ionic  network,  the  network  provides  additional  strength  (slightly  higher 
than  7%  ionic  content)  until  s  0.3  and  then,  stress  drops  rapidly  to  a  level  seen  for  the 
standard  network. 
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Fig.  4  Views  of  the  network  during  deformation  in  the  tensile  mode  at  £  =  0.0,  0.5,  1.0,  1.5,  and  2.0 


Fig.  5  Number  of  bonds  (V^) 
(upper panel)  and  stress  (a) 
versus  strain  (s)  upon  the  shear 
deformation  at  strain  rate  of 
0.6  X  10“^  1/r  for  the  standard 
network  (line)  and  the  ionic 
network  (circles) 


£ 


Snapshots  of  the  network  under  tensile  deformation  at  five  consecutive  strains  are  pre¬ 
sented  in  Fig.  4  to  aid  visualizing  the  deformation  in  tensile  mode.  At  £  =  0.5,  there  are 
voids  formed  in  the  network.  They  actually  emerge  as  the  strain  passes  through  the  yield 
point.  Their  formation  does  not  require  breaking  of  bonds  since  they  are  already  present 
but  are  compacted.  Moreover,  their  presence  indicates  that  the  cross-linking  is  not  entirely 
homogenous:  there  are  regions  of  highly  dense  cross-linking  and  those  of  lesser  amount  of 
cross-links.  They  basically  act  as  defects  that  initiate  cracks.  Ai  s  =  1.0,  the  voids  grow 
further  and  the  surrounding  areas  approach  their  limits  of  the  stress  they  can  carry  without 
breaking  of  bonds.  Later,  at  £  =  1.5,  the  stretching  strands  that  bear  the  load  become  more 
visible.  At  this  stage,  there  is  a  considerable  degree  of  bond  breaking  and  the  network  starts 
to  break  into  two  pieces.  Finally,  dXs  =  2.0,  the  network  has  nearly  failed  with  only  a  single 
strand  holding  the  newly  created  interfaces. 

The  response  to  the  deformation  in  the  shear  mode  is  presented  in  Fig.  5.  The  initial 
rapid  rise  of  the  stress  and  the  subsequent  slowing  down  at  £  =  0. 1  is  due  to  the  yield  be¬ 
havior.  The  magnitude  of  the  yield  stress  for  the  standard  network  is  0.01uq/(P  which  is  in 


9 


^  Springer 


214 


Mech  Time-Depend  Mater  (2008)  12:  205-220 


Fig.  6  Views  of  the  network  through  x-z  plane  during  the  shear  deformation  at  strain  values  of  (a)  0.0, 
(b)  0.7,  (c)  1.3,  (d)  2.0,  (e)  2.6,  and  (f)  3.3.  The  beads  are  colored  in  stripes  to  help  visualize  the  shear 
deformation 


agreement  with  O.OSwo/^^  found  in  the  earlier  study  (Stevens  2001).  The  yield  stress  in  this 
deformation  mode  stems  from  the  frictional  resistance  between  two  sliding  surfaces  interact¬ 
ing  via  LJ  beads.  After  the  yield  point,  stress  continues  to  rise  first,  at  a  lower  rate,  and  later 
(g  >  ~1 .2)  at  a  higher  rate  until  it  reaches  the  maximum  point  (i.e.  where  the  ultimate  failure 
occurs).  Note  that  the  bond  breaking  continues  even  after  the  maximum  point,  because  the 
failure  results  in  a  rough  interface.  This  interface  is  smoothened  while  being  sheared,  thus 
bond  breaking  continues  until  two  sufficiently  smooth  surfaces  are  attained.  Subsequently, 
the  shearing  would  constitute  shearing  of  two  flat  surfaces  interacting  via  LJ  interactions 
similar  to  the  case  at  the  yield  point.  Therefore,  if  simulation  was  pursued  further  the  stress 
would  eventually  level  off  at  about  the  yield  stress. 

In  the  earlier  studies  (Stevens  2001;  Tsige  et  al.  2004)  the  maximum  point  was  followed 
by  a  sudden  drop  in  stress.  This  drop  is  more  gradual  in  our  study  because  of  the  lower  degree 
of  cross-linking  since  densely  cross-linked  networks  fail  more  abruptly.  The  comparison  of 
the  networks  shows  that  the  stress  increases  in  the  ionic  network  more  rapidly  than  the 
standard  network  and  maintains  its  elevated  level  until  the  point  of  the  maximum  stress. 
Similar  to  what  was  found  in  the  deformation  in  the  tensile  mode,  the  incorporation  of  ionic 
pairs  raises  the  shear  stress  thus  generating  a  tougher  network.  However,  after  the  maximum 
stress  is  reached  both  curves  become  nearly  overlapping.  We  discuss  this  behavior  in  the 
following  section. 

The  deformation  the  in  the  shear  mode  is  illustrated  in  Fig.  6(a)  through  (f)  using 
six  snapshots  taken  at  ascending  strain  values.  Ai  s  =  0.7,  the  stripes  are  stretched  to  a 
sigmoidal-like  shape.  As  the  network  is  stretched  further  at  s  =  1.3,  2.0  and  2.6  they  get 
straighter  and  thinner  showing  that  the  strands  are  extended  at  or  beyond  their  limits  before 


^  Springer 


10 


Mech  Time-Depend  Mater  (2008)  12:  205-220 


215 


breaking.  Finally,  at  s  =  3.3,  close  observation  reveals  that  the  location  of  interfacial  failure 
is  near  the  middle  of  the  lower  half  of  Fig.  6(f). 


4  Discussion 

Our  generic  coarse-grained  model  for  the  networks  reproduced  the  general  features  in  a 
stress-strain  curve  such  as  the  yield  behavior,  the  leveling  of  the  stress  after  the  yield  point, 
the  subsequent  strain  hardening,  and  the  failure.  Additionally,  we  presented  some  agree¬ 
ments  to  the  universal  empirical  models  for  thermosets  (i.e.  how  Gy  changes  with  T  and 
how  it  relates  to  E).  Even  though,  the  quantities  we  calculate  are  in  the  LJ  units,  we  attempt 
to  relate  them  to  some  empirical  results.  The  issue  is  to  relate  the  units  of  stress  from  the 
LJ  units  in  our  simulations  to  the  real  units  in  experiments.  The  unit  stress,  uold^,  was  es¬ 
timated  to  be  equivalent  to  43  MPa  by  calculating  the  stress  required  to  separate  a  model 
self-assembled  monolayer  (composed  of  beads  with  the  same  interaction  parameters  as  in 
this  work)  from  a  rigid  surface  to  its  experimental  equivalence  (Stevens  2001).  We  use  this 
approximate  value  in  order  to  put  the  predictions  of  this  model  into  experimental  perspective 
for  a  fully  cross-linked  network  with  no  ions. 

The  estimated  values  for  Oy  and  E  from  the  simulations  were  compared  to  the  experi¬ 
mental  data  (Fig.  7).  Gy  has  a  linear  dependency  onTg  —  T  for  both  the  experimental  and  the 
simulation  data,  as  mentioned  previously.  However,  the  model  predicts  a  lower  slope  than 
the  experiments  and  the  model  fails  to  predict  that  Gy  should  vanish  as  T  approaches  to  Tg, 
as  suggested  by  the  empirical  data.  Nonetheless,  the  model  appears  to  be  providing  a  good 
quantitative  agreement  (with  the  assumption  that  =  43  MPa)  at  about  Tg  —  T  =  125  K 
(i.e.  at  r  =  0.3 r*  in  LJ  units).  Coincidentally,  this  is  the  temperature  when  the  simulations 
of  the  tensile  and  the  shear  deformations  were  performed.  The  estimates  of  the  modulus 


Fig.  7  (a)  Yield  stress  versus 
Tg  —  T .  Estimated  values  for  Oy 
from  the  simulations  {closed 
circle)  were  compared  with 
experiments  for  different  epoxy 
thermosets:  DGEBA/DDM 
{diamond),  DGEBA/BA 
{square),  and  DGEBA/DAO 
{triangles)  (Cook  et  al.  1998).  Tg 
for  the  model  network  of  0.5  T* 
was  assumed  to  be  equivalent  to 
400  K  in  real  units,  (b)  The 
Young’s  modulus  versus  T /Tg. 
Estimated  values  for  E  from  the 
simulations  {closed  circle)  were 
compared  with  experiments  for 
different  epoxy  thermosets: 

RTM6  {square),  Araldite-F 
{diamond)  (Deng  et  al.  2007), 
DGEBA  {triangle  up),  TGAP 
{triangle  left),  and  TGDDM 
{triangle  down)  (Becker  et  al. 
2003).  The  slopes  for  the  fits  are 
2.0  and  2.3  for  the  simulation  and 
the  experimental  data, 
respectively 
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Fig.  8  Views  (through  x-z  plane)  of  pairs  of  beads  forming  bonds  that  eventually  brake  after  applying 
deformation.  The  upper  half  shows  views  for  the  control  network  after  the  deformation  in  (a)  tensile  and 
(b)  shear  mode.  Similarly,  the  lower  half  shows  views  for  the  ionic  network  after  the  deformation  in  (c)  tensile 
and  (d)  shear  mode.  The  lines  highlight  approximate  locations  of  the  ultimate  interfacial  failure 


from  the  simulations  show  a  quite  similar  dependency  on  T /Tg  and  they  are  about  1  GPa 
lower  compared  to  the  experimental  data.  Note  that  our  model  does  not  take  into  account  the 
chemical  structure  of  the  repeating  units  in  epoxy  resins  and  that  the  experimental  control 
parameters  such  as  the  deformation  rate  and  the  sample  size  are  several  orders  of  magnitude 
higher  than  the  simulations.  Therefore,  this  minimal  model  shows  a  remarkable  agreement 
with  experiments  considering  the  simplicity  of  the  model  and  the  differences  in  such  exper¬ 
imental  parameters. 

The  number  of  bonds  that  broke  after  deformation  is  an  important  quantity  as  it  relates 
to  the  failure  mechanism.  We  found  consistently  in  both  the  tensile  and  shear  modes  that 
the  ionic  network  requires  a  larger  number  of  bonds  to  be  broken  before  the  ultimate  failure 
occurs.  In  order  to  understand  the  origins  of  this  behavior,  we  first  determined  the  positions 
of  the  failing  bonds.  These  bonds,  as  identified  in  Fig.  8,  indicate  the  location  of  the  in¬ 
terfacial  failure.  For  example,  the  clustering  of  these  broken  bonds  in  the  lower  portion  of 
the  Fig.  8(a)  points  to  the  location  of  failure  as  confirmed  by  Fig.  4  for  the  deformation  of 
the  standard  network  in  the  tensile  mode.  Figure  8(c)  shows  that,  in  the  case  of  the  ionic 
network,  these  bonds  are  present  not  only  in  the  location  of  the  failure  but  also  in  the  sur¬ 
rounding  area  by  a  larger  quantity  than  the  control  network  (Fig.  8(a)).  This  result  is  even 
clearer  in  shear  mode  (Fig.  8(b)  vs.  (d)). 

We  infer  from  these  observations  that  the  integration  of  ionic  pairs  brings  about  a  network 
that  can  distribute  the  stress  more  uniformly,  which  in  turn  results  in  improved  mechanical 
properties.  In  order  to  gain  insight  into  how  the  stress  is  distributed  more  uniformly  we 
looked  at  the  growth  of  voids  that  form  with  tensile  deformation.  It  is  an  important  phe- 
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Fig.  9  Snapshots  of  the  voids  as  they  develop  during  the  tensile  deformation  in  the  standard  {upper  panel) 
and  the  ionic  {lower  panel)  networks.  Strains  at  which  snapshots  were  generated  are  {from  left  to  right):  0.1, 
0.15,  0.25,0.5,  1.0  and  1.5 


nomenon  because  the  strands  surrounding  voids  are  where  the  stresses  concentrate.  The 
snapshots  of  voids  are  shown  in  Fig.  9  at  varying  strains.  At  the  yield  point  {s  =  0.1),  the 
voids  are  visible  in  their  early  form.  They  appear  to  be  more  randomly  distributed  initially 
in  the  standard  than  the  ionic  network,  while  at  strains  £  >  1.0  they  are  equally  randomly 
distributed.  The  variation  in  their  distribution  is  not  expected  to  affect  the  mechanical  prop¬ 
erties  significantly,  because  the  difference  in  randomness  is  seen  only  until  bond  breaking 
starts.  However,  it  is  notable  to  observe  that  the  sizes  of  the  voids  are  smaller  in  the  ionic 
network  compared  to  the  standard  network  at  any  strain.  The  smaller  void  sizes  the  larger 
the  area  separating  the  material  from  the  empty  space.  This  means  that  the  stress  concentra¬ 
tion  is  not  as  severe  and  the  load  can  be  more  uniformly  distributed  throughout  the  network. 
This  observation  supports  the  analysis  of  the  distribution  of  broken  bonds  in  showing  how 
the  mechanical  properties  are  improved. 

We  had  observed  earlier  that  the  shear  stress  is  enhanced  up  to  the  point  of  the  interfacial 
failure  (Fig.  5).  After  this  point,  stresses  became  indistinguishable  between  the  two  types  of 
networks.  Unlike  this  observation,  when  a  lower  deformation  rate  (by  an  order  of  magnitude) 
is  applied  (see  Fig.  10),  the  ionic  network  continues  to  bear  a  higher  load  even  after  the  point 
of  the  maximum  stress.  Knowing  that  the  two  experiments  differ  only  by  the  deformation 
rate  and  that  the  networks  differ  only  by  the  presence  of  ionic  pairs,  this  result  indicates  that 
it  is  a  consequence  of  the  rate  dependence  of  ionic  pairs.  A  plausible  explanation  is  that  the 
ionic  pairs  that  deform  as  a  result  of  shearing  are  able  to  re-form  only  when  the  system  is 
deformed  with  a  low  strain  rate  so  that  the  strands  have  enough  time  to  relax.  We  infer  also 
that  this  time  dependency  refers  to  the  presence  of  a  self-healing  mechanism  that  restores 
the  some  of  the  strength  of  the  network. 

The  ionic  interactions  were  introduced  to  the  network  in  randomly  distributed  pairs  as 
explained  above.  This  approach  has  been  followed  in  order,  besides  its  practicality,  to  eval¬ 
uate  the  concept  of  ionic  interactions  when  they  are  incorporated  in  the  network  ideally: 
each  charged  particle  interacting  with  its  counterpart.  However,  if  the  ionic  interactions 
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Fig.  10  Stress  (cr)  versus  strain 
(g)  upon  the  shear  deformation  at 
strain  rate  of  0.6  x  10“"^  1  /r  for 
the  standard  network  (line)  and 
the  ionic  network  {circles) 


were  present  during  the  cross-linking  stage,  the  resulting  network  might  have  had  a  dif¬ 
ferent  structure.  While  one  would  expect  to  observe  oppositely  charged  particles  interacting 
thus  being  in  close  proximity,  they  could  be  in  the  form  of  small  clusters  rather  than  pairs. 
Therefore,  the  results  of  this  work  should  be  interpreted  with  this  fact  in  consideration. 


5  Concluding  remarks 

We  presented  results  for  a  coarse-grained  model  of  cross-linked  networks  with  and  without 
ionic  bonds  incorporated.  The  MD  simulations  of  the  deformation  of  the  networks  in  tensile 
and  shear  mode  showed  that  the  model  has  captured  the  fundamental  mechanical  features 
of  thermosets  such  as  the  yield  behavior,  the  leveling  of  the  stress  after  the  yield  point,  the 
subsequent  strain  hardening,  and  the  failure.  Moreover,  we  presented  some  agreements  to 
the  universal  empirical  models  for  thermosets  at  both  qualitative  and  quantitative  level. 

One  is  bound  by  the  length  scale  that  can  be  practically  modeled,  as  an  inherent  limitation 
of  molecular  level  simulations.  This  constraint  reflected  itself  for  the  current  system  in  two 
ways.  First,  how  the  stress  is  transferred  from  the  solid  layers  to  the  network  should  be  a 
function  of  the  distance  between  the  upper  and  lower  layers.  The  interfacial  failures  that 
occur  in  the  vicinity  of  the  layers,  rather  than  in  the  middle  section  in  all  four  cases  may 
be  a  consequence  of  the  stress  build  up  in  this  region.  Secondly,  the  network  can  be  only 
as  heterogeneous  as  the  system  size.  In  fact,  the  main  reason  why  the  system  fails  only 
after  very  high  strains  {s  >  2)  is  that  the  crack  initiation  can  occur  only  after  strands  are 
extended  significantly.  In  a  real  material  there  would  be  a  heterogeneous  cross-linking  (i.e. 
broad  distribution  of  degree  of  freedom  of  strands  to  stretch).  When  the  deformation  is 
applied,  the  strands  with  a  low  degree  of  freedom  to  stretch  will  be  stressed  the  most  and 
break  first.  These  regions  where  such  constrained  strands  concentrate  will  act  as  defects 
that  initiate  cracks  before  the  rest  could  collectively  stretch.  For  this  reason,  the  smaller  the 
samples  size  the  lower  the  probability  of  having  defects  in  the  material  that  results  in  higher 
stress  at  break  (Pascault  et  al.  2002).  From  this  point  of  view,  the  current  model  appears  to 
represent  an  ideal  network  compared  to  a  real  one.  Another  interpretation  is  that  the  model  is 
capturing  the  deformation  behavior  at  or  near  the  crack  tip  as  suggested  previously  (Stevens 
2001;  Tsige  et  al.  2004;  Tsige  and  Stevens  2004). 

The  primary  outcome  of  this  work  is  the  finding  that  the  integration  of  ionic  pairs  in 
an  ideal  manner  elevated  the  level  of  stress  that  the  network  can  withstand  over  a  wide 
range  of  strain  values  in  both  the  tensile  and  shear  modes.  This  implies  that  the  ionic  pairs 
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create  a  stronger  and  tougher  network  compared  to  the  control  network.  We  identified  two 
plausible  mechanisms  through  which  the  ionic  pairs  became  effective.  They  enable  a  more 
uniform  distribution  of  the  stress  through  the  network  and  they  are  able  to  reform  after  being 
displaced  by  deformation  at  low  strain  rates. 

In  future  work,  we  plan  to  investigate  the  effects  of  the  cross-link  density,  content  of 
ionic  pairs  and  the  interplay  between  them.  The  formation  of  the  network  in  the  presence  of 
ionic  reactants  will  also  be  carried  out  to  determine  whether  it  leads  to  a  different  network 
structure. 
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